package fem2;

import inf.jlinalg.IVector;
import inf.jlinalg.SolveFailedException;
import inf.text.ArrayFormat;

import java.io.IOException;
import java.util.Collection;
import java.util.LinkedList;
import java.util.List;

import junit.framework.Assert;
import math2.FunctionRnToR;
import math2.RampFunctionRToR;

import org.junit.Test;
import org.junit.runner.RunWith;
import org.junit.runners.Parameterized;
import org.junit.runners.Parameterized.Parameters;

import fem2.analysis.Analysis;
import fem2.analysis.LinearStaticAnalysis;
import fem2.element.BBarStructuralElement;
import fem2.element.EASStructuralElement;
import fem2.element.StructuralElement;
import fem2.enu.EASModeType;
import fem2.enu.EchoLevelType;
import fem2.enu.IntegrationType;
import fem2.enu.State;
import fem2.enu.StressType;
import fem2.material.StVenantKirchhoffMM;
import fem2.pre_and_post.GidPostStructural;
import fem2.strategies.Strategy;

@RunWith(Parameterized.class)
public class PatchTestQ4 {

	@Parameters
	public static Collection<Object[]> parameters() {
		List<Object[]> params = new LinkedList<Object[]>();

		params.add(new Object[] { 0 });
		params.add(new Object[] { 1 });
		params.add(new Object[] { 2 });
		params.add(new Object[] { 3 });
		params.add(new Object[] { 4 });
		return params;
	}

	private int option;

	/**
	 * 
	 * @param option
	 * <br>
	 *            0: StructuralElement<br>
	 *            1: BBarStructuralElement<br>
	 *            2: Enhanced assumed strain element with 4 enhanced modes<br>
	 *            3: Enhanced assumed strain element with 5 enhanced modes<br>
	 *            4: StructuralElement using selective reduced integration<br>
	 */
	public PatchTestQ4(int option) {
		this.option = option;
	}

	public Mesh createMesh() {
		/*
		 * mesh
		 */
		Mesh mesh = new Mesh();

		double l = 1;
		double h = 1;
		// double rx = Math.random();
		// double ry = Math.random();
		double rx = 0.5; // to avoid any concave quadrilateral
		double ry = 0.5;// to avoid any concave quadrilateral
		double r1 = Math.random();
		double r2 = Math.random();
		double r3 = Math.random();
		double r4 = Math.random();

		// double rx = 0.12151610627091536;
		// double ry = 0.7570319480193257;
		// double r1 = 0.4648281034532379;
		// double r2 = 0.5391821472246238;
		// double r3 = 0.6296424541730731;
		// double r4 = 0.1232934495808714;

		// Debugger.watch("rx = ", rx);
		// Debugger.watch("ry = ", ry);
		// Debugger.watch("r1 = ", r1);
		// Debugger.watch("r2 = ", r2);
		// Debugger.watch("r3 = ", r3);
		// Debugger.watch("r4 = ", r4);

		Node n1 = mesh.addNode(new Node(0, 0));
		Node n2 = mesh.addNode(new Node(l, 0));
		Node n3 = mesh.addNode(new Node(l, h));
		Node n4 = mesh.addNode(new Node(0, h));
		Node n5 = mesh.addNode(new Node(rx * l, ry * h));
		Node n6 = mesh.addNode(new Node(r1 * l, 0));
		Node n7 = mesh.addNode(new Node(l, r2 * h));
		Node n8 = mesh.addNode(new Node(r3 * l, h));
		Node n9 = mesh.addNode(new Node(0, r4 * h));

		mesh.addFace(new Face(n1, n6, n5, n9));
		mesh.addFace(new Face(n6, n2, n7, n5));
		mesh.addFace(new Face(n5, n7, n3, n8));
		mesh.addFace(new Face(n5, n8, n4, n9));

		mesh.addEdge(new Edge(n2, n7));
		mesh.addEdge(new Edge(n3, n7));

		mesh.setNumberOfFields(2);

		/*
		 * constraint
		 */
		Constraint c1 = new Constraint(false, false);
		Constraint c2 = new Constraint(true, false);
		Constraint c3 = new Constraint(false, true);
		n1.setConstraint(c1);
		n2.setConstraint(c2);
		n6.setConstraint(c2);
		n4.setConstraint(c3);
		n9.setConstraint(c3);

		return mesh;
	}

	public Model createModel() {
		Mesh mesh = createMesh();
		/*
		 * element and model
		 */
		double E = 2;
		double nu = 0.3;
		double t = 1.0;
		State ss = State.PLANE_STRESS;
		MaterialModel mm = new StVenantKirchhoffMM(E, nu, t, 0, ss);

		Model m = new Model(mesh);
		for (int i = 0; i < mesh.countFaces(); i++) {
			Element e;
			Face f = mesh.getFace(i);

			switch (option) {
			case 0:
				e = new StructuralElement(f, mm);
				break;
			case 1:
				e = new BBarStructuralElement(f, mm);
				break;
			case 2:
				e = new EASStructuralElement(f, mm, EASModeType.Q1E4);
				break;
			case 3:
				e = new EASStructuralElement(f, mm, EASModeType.Q1E5);
				break;
			case 4:
				e = new StructuralElement(f, mm);
				((StructuralElement) e)
						.setIntegrationType(IntegrationType.SELECTIVE_REDUCED_INTEGRATION);
				break;
			default:
				throw new Error("invalid option");
			}

			m.addElement(e);
		}

		return m;
	}

	@Test
	public void testSolve1() throws SolveFailedException {
		Model m = createModel();
		double F = 1.0;

		double y = m.getMesh().getNode(6).getPosition()[1];

		Force l1 = new Force(m.getMesh().getNode(1), F * y / 2, 0);
		Force l2 = new Force(m.getMesh().getNode(2), F * (1 - y) / 2, 0);
		Force l3 = new Force(m.getMesh().getNode(6), F / 2, 0);
		m.addLoad(l1);
		m.addLoad(l2);
		m.addLoad(l3);

		Analysis an = new LinearStaticAnalysis(m);
		an.getStrategy().setEchoLevel(EchoLevelType.OUTPUT);
		an.run();

		/*
		 * print the displacement and check with theoretical values
		 */
		double[][] Uh = new double[][] { { 0, 0 }, { 0.5, 0 }, { 0.5, -0.15 }, { 0, -0.15 } };
		for (int i = 0; i < 4; i++) {
			double[] uHatLocal = m.getMesh().getNode(i).getUHatLocal();
			System.out.println(ArrayFormat.format(uHatLocal));
			Assert.assertEquals(Uh[i][0], uHatLocal[0], 1e-10);
			Assert.assertEquals(Uh[i][1], uHatLocal[1], 1e-10);
		}

		for (int i = 0; i < m.countElements(); i++) {
			StructuralElement e = (StructuralElement) m.getElement(i);
			for (int j = 0; j < e.getNumIntegrationPoints(); j++) {
				IVector sigma = e.getStress(j);
				Assert.assertEquals(1.0, sigma.get(0), 1e-10);
				Assert.assertEquals(0.0, sigma.get(1), 1e-10);
				Assert.assertEquals(0.0, sigma.get(2), 1e-10);
			}
		}
	}

	@Test
	public void testSolve2() throws SolveFailedException {
		Model m = createModel();
		double F = 1;
		Mesh mesh = m.getMesh();

		FunctionRnToR ld = new RampFunctionRToR(-1, 1, F, F);
		Load el1 = new DistributedLoad(mesh.getEdge(0), ld, null);
		Load el2 = new DistributedLoad(mesh.getEdge(1), ld, null);
		m.addLoad(el1);
		m.addLoad(el2);

		Analysis an = new LinearStaticAnalysis(m);
		an.run();

		/*
		 * print the displacement and check with theoretical values
		 */
		double[][] Uh = new double[][] { { 0, 0 }, { 0.5, 0 }, { 0.5, -0.15 }, { 0, -0.15 } };
		for (int i = 0; i < 4; i++) {
			double[] uHatLocal = m.getMesh().getNode(i).getUHatLocal();
			System.out.println(ArrayFormat.format(uHatLocal));
			Assert.assertEquals(Uh[i][0], uHatLocal[0], 1e-10);
			Assert.assertEquals(Uh[i][1], uHatLocal[1], 1e-10);
		}

		for (int i = 0; i < m.countElements(); i++) {
			StructuralElement e = (StructuralElement) m.getElement(i);
			for (int j = 0; j < e.getNumIntegrationPoints(); j++) {
				IVector sigma = e.getStress(j);
				Assert.assertEquals(1.0, sigma.get(0), 1e-10);
				Assert.assertEquals(0.0, sigma.get(1), 1e-10);
				Assert.assertEquals(0.0, sigma.get(2), 1e-10);
			}
		}
	}

	@Test
	public void testSolve3() throws SolveFailedException {
		Model m = createModel();
		Mesh mesh = m.getMesh();

		double du = 0.5;
		Constraint c1 = new Constraint(false, false);
		c1.setValue(0, du);
		mesh.getNode(1).setConstraint(c1);
		Constraint c2 = new Constraint(false, true);
		c2.setValue(0, du);
		mesh.getNode(2).setConstraint(c2);
		mesh.getNode(6).setConstraint(c2);

		Analysis an = new LinearStaticAnalysis(m);
		an.run();

		/*
		 * print the displacement and check with theoretical values
		 */
		double[][] Uh = new double[][] { { 0, 0 }, { 0.5, 0 }, { 0.5, -0.15 }, { 0, -0.15 } };
		for (int i = 0; i < 4; i++) {
			double[] uHatLocal = m.getMesh().getNode(i).getUHatLocal();
			System.out.println(ArrayFormat.format(uHatLocal));
			Assert.assertEquals(Uh[i][0], uHatLocal[0], 1e-10);
			Assert.assertEquals(Uh[i][1], uHatLocal[1], 1e-10);
		}

		for (int i = 0; i < m.countElements(); i++) {
			StructuralElement e = (StructuralElement) m.getElement(i);
			for (int j = 0; j < e.getNumIntegrationPoints(); j++) {
				IVector sigma = e.getStress(j);
				Assert.assertEquals(1.0, sigma.get(0), 1e-10);
				Assert.assertEquals(0.0, sigma.get(1), 1e-10);
				Assert.assertEquals(0.0, sigma.get(2), 1e-10);
			}
		}
	}

	public static void main(String[] args) throws SolveFailedException, IOException {
		Model m = new PatchTestQ4(0).createModel();

		Mesh mesh = m.getMesh();
		double du = 0.5;
		Constraint c1 = new Constraint(false, false);
		c1.setValue(0, du);
		mesh.getNode(1).setConstraint(c1);
		Constraint c2 = new Constraint(false, true);
		c2.setValue(0, du);
		mesh.getNode(2).setConstraint(c2);
		mesh.getNode(6).setConstraint(c2);

		// new Viewer(m).setVisible(true);

		Analysis an = new LinearStaticAnalysis(m);
		Strategy s = an.getStrategy();
		s.setEchoLevel(EchoLevelType.OUTPUT, EchoLevelType.AUX1);
		an.getModel().setEchoLevel(EchoLevelType.AUX1);
		an.run();

		String projDir = "/home/hbui/kratos_janosch";
		String projName = "test-gid.gid";

		GidPostStructural post = new GidPostStructural(projDir, projName);
		post.clearAllPostFiles();
		post.writeMesh(m, 0);
		post.writeGaussPointStressResults(m, 0, StressType.CAUCHY_STRESS);
	}
}
